import matplotlib.pyplot as plt
import numpy
import pandas as pd
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as py
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
from ipywidgets import interact, SelectionSlider
pio.renderers.default = "notebook_connected"
py.offline.init_notebook_mode()

Numerical Analysis of 2D Diffusion-Reaction Equations Using an Implicit-Explicit Method: Solution and Visualization

• Overview: This is a showcase of a hybrid implicit-explicit method to solve the diffusion-reaction PDE in 2-dimensions. It uses Plotly to provide a stroboscopic and interactive visualization of the solution.

• Connection to data science: The PDE has indirect connections to time-series analysis and forecasting. The diffusion equation can be used to smooth time-series data, removing short-term fluctuations while preserving long-term trends. This is analogous to applying a low-pass filter. The diffusion helps in removing short-term fluctuations, and the reaction term can introduce additional effects such as trends or external influences.

• Why use a hybrid method?: Explicit methods are generally faster than implicit methods, except when the timestep needs to be very small to ensure stability. The diffusion term involves second spatial derivatives, which tend to amplify numerical errors in explicit methods like Forward Euler. The stability criterion for the diffusion term is $\Delta t=\frac{\Delta x^2}{4D}$, where $D$ is the diffusion coefficient. This requires small time steps, making the method not time-efficient. The second derivative term is solved implicitly, allowing for larger time steps without compromising stability, while the reaction term is handled explicitly, reducing the method's time complexity. The hybrid implicit-explicit method is faster than using either a purely implicit or purely explicit method

The Implicit-Explicit Method:

The PDE is

\begin{aligned} u_t &= \sigma D_1 (\partial_{xx}+\partial_{yy})u + f(u, v) \\ v_t &= \sigma D_2 (\partial_{xx}+\partial_{yy})v + g(u, v) \end{aligned}

with the reaction terms

\begin{aligned} f(u,v) &= \alpha u (1 - \tau_1 v^2) + v (1 - \tau_2 u) \\ g(u,v) &= \beta v + \alpha \tau_1 u v^2 + u (\gamma + \tau_2 v). \end{aligned}

The following numerical method solves the PDE over the square domain $\{(x,y):-1\leq x,y\leq1\}$ with periodic boundary conditions.

Numerical method:

\begin{aligned} U^\ast &= U^n + \Delta t \sigma D_1 (\partial_{xx}+\partial_{yy}) U^\ast \\ V^\ast &= V^n + \Delta t \sigma D_2 (\partial_{xx}+\partial_{yy}) V^\ast \\ U^{n+1} &= U^\ast + \Delta t f(U^\ast, V^\ast) \\ V^{n+1} &= V^\ast + \Delta t g(U^\ast, V^\ast). \end{aligned}

The following function discretizes the spacial derivative using a five-point stencil in 2d and returns a sparse matrix representation.

def laplacian_discretization(m):
    """Constructs a sparse matrix that discretizes the 2d Laplacian
    Uses a five-point stencil and periodic boundary conditions.
    """
    delta_x = 2.0 / (m + 1)
    
    # Primary discretization
    e = numpy.ones(m)
    T = sparse.spdiags([e, -4.0 * e, e], [-1, 0, 1], m, m)
    S = sparse.spdiags([e, e], [-1, 1], m, m)
    I = sparse.eye(m)
    A = sparse.kron(I, T) + sparse.kron(S, I)
    
    # Construct periodic BCs
    e = numpy.ones(m**2)
    A_periodic = sparse.spdiags([e, e],[m - m**2, m**2 - m], m**2, m**2).tolil()
    
    # Left & right BCs:
    for i in range(m):
        A_periodic[i * m, (i + 1) * m - 1] = 1.0
        A_periodic[(i + 1) * m - 1, i * m] = 1.0
    
    # Combine two matrices
    A = A + A_periodic
    A /= delta_x**2
    A = A.todia()
    
    return A

def f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
    return alpha * U * (1.0 - tau_1 * V**2) + V * (1.0 - tau_2 * U)

def g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
    return beta * V * (1.0 + alpha * tau_1 / beta * U * V) + U * (gamma + tau_2 * V)

def forward_euler_step(U, V, delta_t, A, sigma, f, g, D1=0.5, D2=1.0):
    """Take a single forward Euler step on the reaction-diffusion equation"""
    
    U_new = U + delta_t * (sigma * D1 * A * U + f(U, V))
    V_new = V + delta_t * (sigma * D2 * A * V + g(U, V))
    
    return U_new, V_new


def backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2):
    """Take a single backward Euler step on the reaction-diffusion equation"""
    U = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_1 * A), U)
    V = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_2 * A), V)
    return U, V

def imex_solver(sigma, tau_1, tau_2, alpha, beta, gamma, D_1, D_2):
    # Alias reaction functions with the above parameters
    f = lambda U, V: f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)
    g = lambda U, V: g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)

    # Set up grid
    m = 150
    delta_x = 2.0 / m
    x = numpy.linspace(-1.0, 1.0, m)
    y = numpy.linspace(-1.0, 1.0, m)
    Y, X = numpy.meshgrid(y, x)

    # Initial data
    U = numpy.random.randn(m, m) / 2.0
    V = numpy.random.randn(m, m) / 2.0
    #fig = plt.figure()
    #axes = fig.add_subplot(1, 1, 1, aspect='equal')
    #plot = axes.pcolor(x, y, U, cmap=plt.get_cmap("viridis"))
    #fig.colorbar(plot)

    # Setup spatial discretization
    U = U.reshape(-1)
    V = V.reshape(-1)
    A = laplacian_discretization(m)

    # Time
    t = 0.0
    t_final = 30.0
    delta_t = delta_x / (10.0 * sigma)
    num_steps = int(numpy.round(t_final / delta_t))

    data = []
    data.append([x,y,U, t])
    
    # Evolve in time
    next_output_time = 0.0
    for j in range(num_steps):
        U, V = backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2)
        U, V = forward_euler_step(U, V, delta_t, A, sigma, f, g)   
        t += delta_t

        if t >= next_output_time:
            next_output_time += 5.0
            U_output = U.reshape((m, m))

            #fig = plt.figure()
            #axes = fig.add_subplot(1, 1, 1, aspect='equal')
            #plot = axes.pcolor(x, y, U_output, cmap=plt.get_cmap("viridis"))
            #fig.colorbar(plot)
            #axes.set_title("t = %s" % t)
            data.append([x,y,U_output, t])

    #plt.show()
    return data

# Parameters
data = imex_solver(sigma=0.0021, tau_1=3.5, tau_2=0, alpha=0.899, beta=-0.91, gamma=-0.899, D_1=0.5, D_2=1.0)
print('checkpoint')

checkpoint

# define time_intervals
time_intervals = []
for i in range(len(data)):
    time_intervals.append(data[i][3])


def extract_data(time):
    j = time_intervals.index(time)
    z_data = data[j][2]
 
    return z_data

def update_plot(time):
    
    z_data = extract_data(time)
    
    z_data = z_data.reshape((150, 150))
    

    fig = go.Figure(data=[go.Surface(z=z_data)])

    
    fig.update_traces(contours_z=dict(show=True, usecolormap=True,
                                      highlightcolor="limegreen", project_z=True),)
    fig.update_layout(title='Topographical 3-D Surface Plot of the Solution with Contours at t={}'.format(time), autosize=False,
                      scene_camera_eye=dict(x=1.7, y=1.3, z=0.64),
                      width=1000, height=700,
                      margin=dict(l=85, r=50, b=65, t=90),
    )

    fig.show()

for i in time_intervals:
    update_plot(i)
    
# These widgets are currently configured for .ipynb environments
# time_slider = SelectionSlider(options=time_intervals, value=time_intervals[5], description='Time')
# interact(update_plot, time=time_slider);

References

[1] Mandli, K. T. (n.d.). Reaction-Diffusion Demo. In Numerical Methods for Partial Differential Equations [Jupyter Notebook]. GitHub. Retrieved July 27, 2024, from https://github.com/mandli/numerical-methods-pdes/blob/master/reaction-diffusion_demo.ipynb

[2] Ketcheson, D. (n.d.). Finite Difference Course. GitHub. Retrieved July 27, 2024, from https://github.com/ketch/finite-difference-course

[3] LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady State and Time Dependent Problems. Society for Industrial and Applied Mathematics (SIAM).